Understanding QMC algorithms

I started working on QMCPACK earlier this year.  As usual when diving into an unfamiliar code base, one of the first challenges is to understand how the code works.   Scientific codes can take more work to comprehend because they implement complex algorithms, and one must understand the underlying algorithms in addition to the code structure.  Reverse-engineering the algorithm from the code is usually challenging.  Hopefully there's an associated paper that describes the algorithm.  However, publications give a high-level view of the algorithm and the mathematics, but there's usually a gap between the paper and the details of the implementation.

Once the code and algorithm are understood, the next question is how to best test the code.   Ideally the original implementer should leave behind code and documentation for this purpose, but usually what we have is the final code and a paper.    What other kinds of artifacts (code, documentation, data, test cases,  etc.) should an implementer produce to help with this?

In an attempt to address these issues, I created a repository with various Python scripts and Jupyter (IPython) notebooks covering different algorithms used in a Quantum Monte Carlo code.   The repository is here: https://github.com/markdewing/qmc_algorithms

Some more benefits and goals of these scripts include

Ease of comprehension

Notebooks show the derivations and demonstrate the algorithm working in a simpler context.  These aim to be a bridge between the paper or text description of an algorithm and the final production implementation.

Symbolic Representation

Represent the various equations in a more symbolic form.  This allows mathematics-aware transformations, such as taking derivatives.   See the Jastrow factor in Wavefunctions/Pade_Jastrow.ipynb for an example.   Also the wavefunctions and expressions for the energy in Variational/Variational_Hydrogen.ipynb and Variational/Variational_Helium.ipynb.

Testing

These scripts can be used for testing in several ways.  One is to evaluate the expressions at a single point to check values in unit tests.

Eventually the symbolic or script representation should be compared against the code for a wider range of values and parameters.  This is similar to property-based testing in software engineering, but there are a couple of issues to address.  The first is infrastructure support.  Most property-based testing frameworks focus on integer parameters, like array lengths, and not so much on floating point values.  A second issue is cross-language coding.  In our case, the script representation is in Python and the code to test is in C++. 

Also for testing, simpler or slower algorithms can be used to verify more complex algorithms.  For example, the repository contains a python script for Ewald summation (LongRange/ewald_sum.py) that uses the standard break-up between short and long range pieces.  QMCPACK optimizes this split between the short and long range pieces, but the result should still be the same, up to some error tolerance.

Code generation

Ultimately, code generation can be used to connect the symbolic form to the executable code. This will a require a maintainable mechanism for doing the code generation (it's easy to start doing code generation, but there can be problems with scaling it and keeping it neat, clean, and understandable. ) The previous step - using the symbolic forms for testing - seems like an easier intermediate path, and is better suited for testing existing code.

Summary

This is an experiment in creating artifacts that enhance comprehension and testing of scientific programs.  Take a look at the repository ( https://github.com/markdewing/qmc_algorithms ) and see if you find it useful.

Reduced-order models in QMC

After the last post about multiscale modeling, let's explore how some of the methods and ideas might apply to improving QMC.

Possible sources:
- Reduced order models.  Usually applied to time-dependent problems.
- Dimensionality reduction.  From statistics and machine learning
- Low-rank matrix approximation.
- Compressed sensing - determine a sparse sampling before performing the sampling.  See this article for an application to wavefunctions.

Some of these approaches are applications to different domains but use similar underlying mathematics (e.g., apply Singular Value Decomposition (SVD) to an appropriate matrix).

The "order" in "reduced-order models" refers to the complexity of the model.  Some possibilities for QMC on how the complexity might be varied include: number of electrons, number of sample points in the integral, quality of the wavefunction, and the precision of the arithmetic.

Let's look at each of these in a little more detail:

1. Number of electrons
   A. Integrate out some electrons - yields pseudopotential
   B. Integrate out all electrons - yields intermolecular potentials between atoms/molecules

 To ultimately reduce the complexity of the model, most or all of the electrons need to be eliminated.  This is also the natural decomposition at an initial look - generate potential from electronic structure calculation.    However, this may be too large of jump of scales in some cases, and it may be useful to find intermediate models that still include some or all of the electrons.

2. Number of sample points in the integral
    A. Subset of points from an existing set of samples.  Not so useful for that particular run, but could be useful for a nearby system. Possible ways a system could be 'nearby':
         - Changing parameters for VMC optimization
         - Core electrons for psuedopotentials
         - Weak intermolecular interactions
   B. Generate new points faster or more accurately
        Quasi Monte Carlo methods (the other QMC) or other low discrepancy sequences.

Side note on pseudopotentials.  They serve two purposes:
   a. Reduce the number of particles that need to be handled
   b. Reduce the variance introduced by sampling the core electrons

  I suspect that (b) is more important than (a), at least for small problems.  Pseudopotentials complicate various QMC algorithms, so improved sampling that reduces the variance could be an effective alternative.

3. Quality of wavefunction
  In general we would like to have a systematic way of controlling the quality (computational complexity vs. error) of the wavefunction.   One use for less-accurate approximations to the wavefunction might be as a guide for faster generation of new sample points.

4. Precision of arithmetic
  The full model is really the position of all the electrons in space.   We might consider restricting the electrons to lattice points or reducing the precision in places.   For examples in machine learning, see here and here.